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ROBUST PRECONDITIONERS FOR INCOMPRESSIBLE MHD MODELS 


YICONG MA, KAIBO HU, XIAOZHE HU, AND JINCHAO XU 

Abstract. In this paper, we develop two classes of robust preconditioners for the structure¬ 
preserving discretization of the incompressible magnetohydrodynamics (MHD) system. 
By studying the well-posedness of the discrete system, we design block preconditioners 
for them and carry out rigorous analysis on their performance. We prove that such pre¬ 
conditioners are robust with respect to most physical and discretization parameters. In 
our proof, we improve the existing estimates of the block triangular preconditioners for 
saddle point problems by removing the scaling parameters, which are usually difficult 
to choose in practice. This new technique is not only applicable to the MHD system, 
but also to other problems. Moreover, we prove that Krylov iterative methods with our 
preconditioners preserve the divergence-free condition exactly, which complements the 
structure-preserving discretization. Another feature is that we can directly generalize this 
technique to other discretizations of the MHD system. We also present preliminary nu¬ 
merical results to support the theoretical results and demonstrate the robustness of the 
proposed preconditioners. 


1. Introduction 


The incompressible Magnetohydrodynamics (MHD) system models the interactions 
between electromagnetic fields and conducting fluids. It consists of the incompress¬ 
ible Navier-Stokes equation for the fluids and the (reduced) Maxwell's equation for the 
electro-magnetic fields. MHD systems of different scales are used in different fields, 
such as astrophysics, engineering related to liquid metal, controlled thermonuclear fu¬ 
sion. There is a vast literature on the study of various aspects of MHD systems. In this 
work, we concentrate on the following incompressible MHD system in both 2D and 3D. 
We assume that Cl c 1R 2 or Q c 1R 3 is a simply connected bounded domain with a Lipchitz 
boundary. In 3D case, the model is 


(1.1) 

du , 1 

-, {U -W)u-- 

(1.2) 

?)B „ 

—- + Vx£ = 0, 
at 

(1.3) 

j ~wC x '‘' lB ■ 

(1.4) 

a r (E + ux B) = j, 

(1.5) 

V • it = 0. 


Here, u is the velocity of fluid, p is the pressure, B is the magnetic field, j is the volume 
current density, and E is the electric field. The physical parameters are the fluid Reynolds 
number Re, the magnetic Reynolds number Rm, the coupling number s, the relative elec¬ 
trical conductivity c r/ and the relative magnetic permeability p r . The initial conditions for 


Key words and phrases, incompressible MHD, robust preconditioners, field-of-values analysis. 

1 



2 


YICONG MA, KAIBO HU, XIAOZHE HU, AND JINCHAO XU 


the fluid velocity, and the magnetic field are 

u(x, 0) = uq(x), B(x,0) = Bq(x), Vx E Cl. 
And the boundary conditions are 

u = 0, n ■ B = 0, n x E - 0, \/x G dCl, t > 0. 


The primary unknown physical variables in the model are u, p and B. These quanti¬ 
ties, once known, uniquely determine E and j. To discretize this system, we follow the 
scheme proposed in our previous work [16], which ensures B is divergence-free exactly 
on the discrete level (Such property is referred to as structure-preserving in the following 
sections). Therefore, we solve u, p, B, and E simultaneously 
In 2D case, the MHD model (1.1)-(1.5) becomes: 


( 1 . 6 ) 

(1.7) 

( 1 . 8 ) 

(1.9) 

( 1 . 10 ) 


du , 1 . 

— + (u ■ v)u — — Aw — sjxB+vp-f, 

dB 

—^r — + curl£ = 0, 
at 


1 -i „ 

] ~R^i mi}lr B = 0 ' 
cr r {E + u x B) = j, 


V • u = 0. 


Here, rotw = —-—1 for any vector u = (u-i , u 2 ) , and curl m = —, — — ) for any scalar 

dx dy \dy dx J 

u. 

Note that we can directly apply the analysis of 3D model to 2D case because we can 
write velocity u = (ui,M 2 ) T e E 2 as m = (iii,ii2,0) T e 1R 3 , the magnetic field B = (B 1 ,B 2 ) t as 
B = (B 1 ,B2,0) t , and the electric field £ as E - (0,0, £) r G 1R 3 . The cross product and de¬ 
rivative operators are all well-defined by rewriting those 2D variables in the 3D fashion. 
Therefore, we focus on the analysis in 3D case in the rest of the paper. 

For the MHD system, solving the linear systems obtained after linearization is usually 
the most challenging and time-consuming part in the overall simulation, which is due to 
the large-scale, multi-physical, and indefinite properties of the resulting linear systems. 
In order to improve the efficiency of the numerical simulations, there have been a lot of 
studies on the development of efficient solvers for various MHD systems. 

Due to the block structure of the resulting linear systems, many block preconditioners 
have been developed in the literature for the MHD system. Shadid and his collaborators 
have developed a series of robust and scalable Newton-Krylov solvers for the MHD sys¬ 
tem [30, 31, 9, 23]. In [31], they propose a robust, efficient, fully-coupled stabilized finite 
element formulation for resistive MHD, which enables both fully-implicit and direct-to- 
steady-state solutions. They investigate the performance of one-level Schwarz method 
and also a new fully coupled algebraic multilevel method in that paper. In [30, 9, 23], they 
explore a class of robust and scalable parallel preconditioners for Newton-Krylov solver 
based on the physical-based approximate block factorization (ABF) technique. They em¬ 
ploy block factorization and approximate the resulting Schur complement recursively 
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based on special techniques, for example, operator commutativity [12]. Numerical exper¬ 
iments and benchmark tests demonstrate the efficiency and scalability of their precon¬ 
ditioners. Chacon and his collaborators also contribute to developing "physics-based" 
block preconditioning strategies for fully-implicit Newton-Krylov solvers for MHD sys¬ 
tem [8, 7, 6, 5]. They use physical-based ABF approach to design preconditioners for 
the linearized MHD system. With algebraic techniques and recursive approximations of 
the Schur complements, they successfully convert the complicated problems into several 
Poisson-like equations and design efficient ABF preconditioners for the linearized MHD 
system. The implementation of ABF preconditioners consists of solving a sequence of 
Poisson-like equations, for which multigrid (MG) methods, especially algebraic multi¬ 
grid (AMG) methods, can be effectively applied. Numerical experiments and benchmark 
tests demonstrate the efficiency and scalability of ABF preconditioners. Moreover, Toth 
et. al [33,17] use a block incomplete LU (ILU) factorization to precondition the MHD sys¬ 
tem. And Badia et. al [1] propose a recursive version block ILU preconditioner recently 

In addition to the block preconditioners, there are also many works on other precondi¬ 
tioning strategies for the MHD system, such as the additive Schwarz methods [25, 4, 21, 
26], Operator splitting method [27, 26]. 

In this paper, we develop robust block preconditioners especially for the linearized 
system arising from the structure-preserving discretizations [16]. We precondition it by 
converting coupled MHD systems into subsystems for which effective preconditioners 
exist. Different from the aforementioned preconditioners that have been studied mostly 
from algebraic point of view, our preconditioners are motivated from the perspective of 
functional and PDE analysis following a framework summarized by Mardal and Winther 
in [20]. In essence, we study the mapping property of the linearized operator between 
appropriate Sobolev spaces equipped with carefully chosen norms. So we can derive 
robust block diagonal preconditioners based on proper norms straightforwardly. Such 
block diagonal preconditioners are often known as norm-equivalent preconditioners and 
use them in combination with minimal residual (MINRES) method. Moreover, we can 
design block triangular preconditioners, and theoretically prove that they are Field-of- 
values- (FOV-) equivalent preconditioners based on the mapping properties [19]. And 
we use them in combination with general minimal residual (GMRES) method. 

In the analysis of FOV-equivalent preconditioners for saddle point problems, we im¬ 
prove the estimates by Loghin and Wathen in [19] by removing scaling parameters in 
front of the diagonal blocks. It is observed that such scaling parameters are difficult to 
choose and unnecessary in practice. By choosing appropriate norms in the analysis, we 
are able to get rid of these scaling parameters, which is consistent with the practical im¬ 
plementations and observations. While this new technique is originally motivated for 
FOV-equivalent preconditioners for the MHD systems, it is expected to be applicable to 
other saddle point type problems. 

One special feature of our work is that we pay special attention to the structure-preserving 
property. We design our preconditioners in such a way that the resulting preconditioned 
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Krylov iterative methods inherit this property even if the linear system is solved inex¬ 
actly This feature makes our preconditioner structure-preserving and suitable for accu¬ 
rate and efficient MHD simulations. 

The efficiency of our preconditioners depends on how the several relevant subsystems 
are solved. The subsystem for velocity field is Poisson-like, for which multigrid methods 
are effective preconditioners. The subsystem for pressure is well-conditioned and hence 
it can be easily solved. For the electric and magnetic fields, we need to solve subsystems 
involving curl curl and grad div operators. We adopt the HX-preconditioner [15], which is 
developed based on the auxiliary space preconditioning [35]. Taking advantage of those 
efficient sub-problem solvers, we develop practical and scalable preconditioners for the 
structure-preserving discretization of the MHD system. 

The rest of the paper is organized as follows. We revisit the structure-preserving fi¬ 
nite element discretization introduced in [16] in §2 and §3. We carry out the analysis 
using different weighted norms to ensure the robustness of preconditioners with respect 
to the physical and discretization parameters. Then we propose and analyze these pre¬ 
conditioners in §4 and discuss their generalizations to other discretization schemes in §5. 
Finally, we present results of numerical experiments in §6 to demonstrate the robustness 
of these new preconditioners. 

2. Magnetohydrodynamics model 

Following [16], we use the following set of notation. First, (■, ■) and || • || denotes L 2 inner 
product and L 2 norm 

{n,v) = / u-vdx, 11 1 /11 = ( / |n| 2 dx) = \/(u,u). 

J n \J n ) 

With a slight abuse of notation, we use L 2 (0) to denote both the scalar and vector L 2 
space. Given a linear operator D, we define 

H(D, Q) = \v G L 2 (0), Dv G L 2 (0)}, 

and 

Hq(D, O) = {v G H(D, O), t^v - 0 on 30} . 

Here, t D is the trace operator defined by 

! v, D = grad, 
v x n, D = curl, 
v ■ n, D = div, 

where n is the outer normal direction of 30. We note that L 2 (0) can be viewed as H(id, O) 
where id denotes the identity operator and we often use the following notation: 

L§(0) = G L 2 (0), J^v = 0 

When D = grad, we often use the notation: 

H 1 ] O) = H(grad, O), H^(O) = H 0 (grad, O). 
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We also use the space U and H 1 with their canonical norms 

Mo, P = ([ M p ) , IMIo,oo=esssup|i7(z)|, |M| h -i= su P 

\‘ /n ' ^eTTl(O) II v II 

As we will see later, it is convenient to introduce the following spaces 

X = H^(n ) 3 x H 0 (div;n) x H 0 (curl;fi), Q = L§(O), 

V = Hg(O) 3 , = H 0 (div; n) and V c = H 0 (curl; D), 

v dfi = H 0 (divO,n) = {c e v d , v ■ c = o}. 


We use W* (W = V, V d or V c ) to denote the dual space of W, and W h the correspond¬ 
ing finite element space of W. Moreover, we assume that both u r and cr r are positive 
continuous functions only depending on x € O, which induce weighted L 2 -norms 


11 33 11 = (' PrX,x ), ||as||f 1 = (}i r 1 x,x). 

For the sake of convenience, we also define a tri-linear form 

1 

d(w, u, v) - - [(to ■ V u, v) — (w ■ V v , u )]. 


Based on these notation, we consider the variational formulation for the incompressible 
MHD system (1.1)-(1.5): Find (u, B,E) e X and p e Q such that for any (v, C,F) € X and 
9 £ 0 / 


( 2 . 1 ) 


( ^-,v ) + d(u, u, v) + k 1 (V • u, V • v) + Vu) 

\ at ) lie 

—s(a r E x B, v) + s(cr r u x B, v x B) — (p, V ■ v) - (/, v ), 
< +(Fr _1 V X E,C) = 0, 

( a r E , F) + (cr r u x B,F) - — — (u7 1 B, V x F) = 0, 

Km 

(V-tM 7 ) = 0 . 


Remark 2.1. Note t/mt f/ie special treatment of the nonlinear convection term in (2.1) is based on 
the following identity, i.e. if V ■ u = 0 and u = 0 on 8 fi, 

(2.2) (u ■ Vm, v) - d(u, u, v), 

This is a classical stabilization technique, c.f [32], 

3. Finite element discretization for the MHD system 

In this section, we revisit the structure-preserving discretization of the incompress¬ 
ible MHD system [16]. For the temporal discretization, we adopt the backward Euler 
scheme. And similar spatial discretization is also applicable to other temporal discretiza¬ 
tion schemes, such as Crank-Nicolson or backward differentiation formula (BDF). We 
first introduce the full discretizations and then revisit the well-posedness of the linearized 
problem with different weighted norms from [16]. 
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3.1. Full discretization scheme. Before going into the details of discretizations, we first 
introduce the finite element spaces for the velocity u, the pressure p, the magnetic field 
B, and the electric field E, respectively 

For the velocity and pressure, we choose a standard stable Stokes pair such that Vu c 
h 0 W, Q h C Lg(fl), and it satisfies the well-known inf-sup condition: 


(3.1) 


inf sup 

qh&Qh Vh eV h 


(V ■ vh,%) 
|V« h || \\cjh II 


>/3>0, 


where the positive constant ( 6 is independent of mesh size h. Many stable Stokes pairs 
are available, such as Taylor-Hood element [2]. For the magnetic field, we use Raviart- 
Thomas elements and denote the finite element space as V)f c H 0 (div; Q). And we employ 
Nedelec edge elements to discretize the electric field and denote the finite element space 
as C Ho (curl; fi). Therefore, we can define the finite element space of X 

X h = V h x V h d x V£. 


Based on the above finite element spaces, we have the following full discretization 
scheme based on the Picard linearization. We discretize the convection term explicitly 
because it leads to a symmetric linear system, which facilitates the analysis of solvers. 
The solvers proposed in this paper also work for other discretizations in [16]. 

Symmetric Picard linearization. Find (u’l Bf t , E" t , p n h ) e X h x Q h such that for any 
(v h ,C h ,F h ,q h ) C Xj } x Q/„ 


(3.2) 


k 1 (u)\-u) 1 i l ,v^j+d(ul 1 ,u n h 1 ,v h ) + k 1 (X-u h ,X-v h ) 

+ Ye <y<,Vv h )-s(jZn -1 >< BF\n) ~ (P^V-« h ) = (fi;,v h ), 
< -*" 1 * (V 1 (B n h - B^- 1 ) , C h ) - a (ji - 1 V X El C h ) = 0, 

s (. jh, n -v F h ) - « [}iV l Bl V x = 0, 

. (V-<,^) = 0. 


where = n r (H ; " + u n h x B^ 1 ), and a = s/Rm. Here, the second equation is multiplied 
by —a and the third equation is multiplied by s. This is because we would like to make 
the resulting linear system symmetric. The above discretization has nice properties, for 
example, energy estimate, structure-preserving, as analyzed in [16]. 


3.2. Well-posedness. Now we discuss about the well-posedness of scheme (3.2), which 
is the foundation of the preconditioners we propose. 

For the sake of simplicity, we rewrite x ( x = u, B, E or p, so is the x mentioned after¬ 
ward in this paragraph) as x~ and x'l as x. We keep the subscript h for the finite element 
spaces in consideration of clarity. Then we can write the symmetric Picard linearization 
as: Find (u, B, E, p) e X^ x Q h such that for any (v, C, F, q) e X h x Q/,, 

k~ 1 (u,v) + Re~ 1 (Vu,'Vv) +W 1 (V • u, V • v) — s(cr r E x B ,v) 

+s(a r u x B~, v x B~) — (p, V ■ v) = (/, v), 

—k~ 1 a.(pp 1 B, C) - Ci{plV x E,C) = -k~ x B~ , C), 
s{a r E, F ) + s(a r u x B ~, F) - tz(}ip 1 B, V x F) - 0, 

(V ■ u,q) - 0, 


(3.3) 
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with f = f + k 1 u - d(u ,u ,v). For £ = (u, B, E) 6 X h , tj = ( v , C, F) e X h , and p, q £ Q/„ 
we define a bilinear form ao(-, ■) on X h x X h and b(-, ■) on X h x Qh by 

a 0 (£, v) = k~ l (u, v) + Re _ 1 (Vit, V«) + fc _1 (V ■ u, V ■ v) — s(cr r E x B~ ,v) 

+ s(u r u x B~,v x £C) — k~ 1 K(pp 1 B, C ) — x F, C) 

+ s(cr r E, F) + s(a r u x F , F) - oc{p~ l B, V x F), 


and 

b {n,q) = (V ■ v,q). 

Therefore, we can write (3.3) as: Find £ e X h and p e Q/, such that 


(3-4) 


ao(€, V) + b (V, V ) = (h, r?), V?/ € X h , 
b (€,q) = (g,?), e Q ft . 


In order to analyze the well-posedness of this problem, we analyze an auxiliary problem 
first. Define a bilinear form 


a(£, v) = a o(£, v) - Oi{p r 1 v ■ B, V ■ C). 

The corresponding auxiliary problem is: Find £ e Xj t and p 6 Q/, such that 

, 3 ^ | «(6 t?) + b (v, V ) = (h, ry), V?/ e X, ; , 

I b (tq) = € Q h . 


As shown in [16], the mixed formulations (3.4) and (3.5) are equivalent if h - ( f,l,r) e 
Vh X [VhT x [*£]* such that (l, C) = (l R , C ), VC e V d for some l R e V d ’°. The well- 
posedness of the mixed formulation (3.4) follows directly from that of the auxiliary prob¬ 
lem (3.5). Therefore, we focus on proving the well-posedness of the auxiliary problem. 

As usual, we use Brezzi's theorem to analyze the well-posedness of the above auxil¬ 
iary problem. As discussed in [20], ensuring that the constants appearing in the inf-sup 
conditions are independent of the physical and discretize parameters is crucial to design 
robust block preconditioners for coupled systems. Therefore, we introduce the weighted 
norms 

\\(v,C,F)\\ 2 x = IMI^ + IIcii^+iifii^, \\q\\ Q = \\ q \\ U2 , 

(3-6) ||(t;, (? ,C,F)||^=||t;||^ i + || (? ||^ 2+ ||C||^ 3 + ||F||^ 4 , 

with 

11^11^= fc~ 1 ||t?|| 2 +Re -1 ||Vr>|| 2 +fc -1 ||V ■ t>|| 2 +s||t> x £C|| 2 r , 
h\\n^k\\q\\ 2 , 

\\c\\ 2 n = fc-^licii^^nv.cii^, 

ll^llw 4 =s||F||^ r +te||VxF|| 2 rl , 

where H, (i = 1, 2, 3 or 4) is a symmetric positive operator (SPD) such that ||a:||^.= ('Hpr,, x). 

Results similar to the following theorem can be found in [16], by introducing more 
carefully chosen weighted norms, this theorem provides better bounds than the previous 


one. 
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Theorem 3.1. Given h e X* t and g e Q* h , the auxiliary problem (3.5) is well-posed if k is 
sufficiently small, i.e. k < k 0/ where 

(3-7) h = l s \W^B-%%. 

That is, for any h e X ; * and g e Q ; *, there exists a unique (£, p) - ( u, B, E, p) e X h x Q h that 
solves the above problem and satisfies: 

\\ u \\ Hl ^B\\ n3+ \\E\\ Hi+ \\p\\ H2 < c (IIMIx’+llgllQ*) ' 

where the constant C does not depend on the mesh size h, the time step size k, and physical parameters Rm, 
s, p, and Vy. 

Proof of Theorem 3.1 is similar to that in [16]. We step by step prove the following 
properties. 

( 1 ) a(-, ■) and b(-, ■) are bounded; 

( 2 ) inf-sup condition holds for bf, ■); 

(3) inf-sup conditions hold for a(-, ■) in the kernel of the operator induced by bf, ■). 
Lemma 3.1. Both a(-,-) and b(-,-) are bounded operators. That is, 

a(tv) < C||£||xIMIx/ 

Hv>q) < C\\r]\\x\\q\\Q, 

where C is a constant independent of the mesh size h, the time step size k, and physical parameters 
Rm, s, p r , and cr r . 

Proof. It is enough to show that every term in a(-, ■) and b(-, ■) is bounded. For the Navier- 
Stokes equation part, 

k~^\(u, v)\ + Re^ 1 \(\7u, Vt>)| +fc _1 |(V ■ u,X7 ■ v)\ < C||m||% 1 ||i;||% 1 . 
with a constant C independent of k and h. For the nonlinear term, 

I s(a r u x B~,E)\ < s||it x B _ || oi .||.E|| (rr < \\u\\ Hi \\E\\ H a , 

and 

\s{(T r U X B~, V X B~) I < S || U X X B \\cr r < || U || ^ || V || ^ . 

Moreover, we have 

k(^r _1 V X E,C )I < a|| V x L;|| ?( -i||C|| ?i -i< \\E\\ Hi \\C\\ Hr 

Note that, 

\s(a r E,F)\<s\\E\\ ar \\F\\ ar < ||25|| W4 ||F|| W4 , 

Ik-Mm'B.Ol < fc- 1 «||S|| ?/ _ 1 ||C7|| ;t -i< ||B|| W 3||C|| W3/ 

and 

|(V ■ v,q)\ < \\v\\ Hl \\q\\ H2 . 


These estimates conclude the lemma. 


□ 
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Lemma 3.2. Ifk<k 0 , defined in (3.7), there exists a constant £ > 0 such that 

inf sup > £ > 0, 

^Q v eX h\\xh\\Q 

and £ zs independent of the mesh size h, the time step size k, and physical parameters Rm, s, y r , 
and cr r . 

Proof. It is well-known that the inf-sup condition of bf, ■) for stable finite element pairs 
holds, i,e. 

inf sup > £ > 0, 

^QvgX IMIiIMIo 

where £ is a constant independent of k. Therefore, for any q C Q, we can choose r/ = 
( v , 0,0) e X such that 

IMIiIMIo 

Note that 

\MhM\h 2 < c IMIiIMM 

and the constant C is independent of k, which completes the proof. 

Lemma 3.3. Ifk<ko, defined in (3.7), there exists a constant f> > 0 such that 

«(£/ v) 


□ 


inf sup -r——fj—I,— > (6 > 0, 

<V£eX°'“ 0 ^ eX o,u iKlIxllMIx 

inf sup ,, — > yS > 0, 

0 ^ eX o,u ||£||xlMlx 

and fS is independent of the mesh size h, the time step size k, and physical parameters Rm, s, y r 
and <r r . Here, 

X°' u = {ue V h , (V ■«,</) = 0, \/q e Q h } . 


Proof. Take v = u, F = E, C = —-(B + kV x E ), then 


ak, 


a (£r r l) = k 1 ||«|| 2+ ^ e 1 ||Vm|| 2 +A: 1 1 | V ■ m|| 2 +s||.E|| 2 + — 1 | V x E\\ 2 ^ 

r 2 Vr 

+ C %-\\ B \\ 2 u -i+M<rrU x B , E) + s||m x B“|| 2 r + ^||V ■ B\\ 2 j. 

Z rr Z rr 


Since 


2 {(JrUxB ,E) <2\\E\\a r \\ux B \\ ar < ^\\E\\ 2 r +2\\u x B 


— 112 
0-r' 


When k < ko, we get s||u x £?~ || 2 r < -fc _1 IMlM Therefore, 

fc _1 ||M|| 2 -s||u x B ~|| 2 r > ^W 1 ||w|| 2 . 

Hence, a(^,r/) > f\\€\\ 2 x - Obviously, we have |M|x< C||£||_x> then the inf-sup condition 
holds. And neither C nor f> depends on k, h, Rm, s, a r or y r . The other inf-sup condition 
can be proved in the same way. 

□ 

Mixed formulations (3.4) and (3.5) are equivalent, and (3.5) is well-posed. As a result, 
we have the following well-posedness for the mixed formulation (3.4). 
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Theorem 3.2. Given h = ( f,l,r ) e Vf x [V d \* x [V£\* such that 

(l,C) = (l R ,C), VC e V d for some l R £ V d '°, 

and g £ Ql, where (•, •) zs f/ze duality pair between V d and [V d \*, the mixed formulation (3.4) is 
well-posed ifk < k 0 . 

Proof. By similar argument in Lemma 1 and Theorem 8 in [16], it is straight-forward to 
reach the conclusion. □ 


4. Robust preconditioner 

In this section, we develop robust preconditioners based on the well-posedness of the 
discrete MHD system. We discuss two types of preconditioners, one is norm-equivalent 
preconditioners, which lead to block diagonal preconditioners, and the other is FOV- 
equivalent preconditioners, which lead to block triangular preconditioners. Here, we 
follow the classification proposed in [19]. 

4.1. Norm-equivalent preconditioner. Following the notation and approaches in [20], 
we consider a general model problem on a Hilbert space H, which is equipped with an 
inner product (■,■)« and an induced norm ll*llw= {x,x)u‘. Find x £ H such that 

(4.1) L{x,y) = (f,y), Vy £ H, 

where L(-, ■) is a symmetric bilinear form, and (■, ■) is the duality pair between H and H*. 
We assume the model problem (4.1) is well-posed and satisfies 

(4.2) \L(x,y)\ < P\\x\\ H \\y\\ H , Vx, y £ H, 

and 

(4.3) inf sup — > a > 0. 

0 ¥x£H 0 -t yeH ||*||?z||y||?z 

We define an operator A: H — >■ H* by 

(Ax,y) = L(x,y), Vx, y £ H. 

Hence, the operator form of (4.1) is 

(4.4) Ax = f. 

Assume that an SPD operator M : H* — >■ H is a preconditioner of this system. Based 
on M, we can define an inner product (x,y) M ~i = (M~ 1 x,y) on H, which induces a norm 
11 * 11 ^- 1 = (x,x) M -i. As 

(MAx,y) M -i = (Ax,y) = L(x,y) = (x,MAy) M - 1 , 

MA : H —>• H is symmetric with respect to (■,-).v( '• Therefore, we can use precondi¬ 
tioned MINRES method [22] to solve (4.4) with M as the preconditioner and (■, -) M ~i as 
the inner product. It is proved [13] that if x m is the zzz-th iteration of MINRES method 
and x is the exact solution, then there exists a constant 5 £ ( 0 , 1 ), only depending on the 
condition number k(M A), such that 

(4.5) {MA{x-x m ),A{x-x m ))^ 2 < 26 m {M A(x - x°), A(x - x 0 )) 1 / 2 . 

Moreover, an estimate leads to 

, _ k(MA) - 1 
“ k(MA) + 1 ' 
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Therefore, in order to estimate the performance of the preconditioner M, we need to 
estimate the condition number k(MA). Mardal and Winther [20] prove that if (4.2) and 
(4.3) hold, and the operator M satisfies 

( 4 . 6 ) ci(x,x) w < ( x,x ) M -i < c 2 (x,x) n , Vx £ H, 

then the condition number satisfies 

1 < k(MA) < —■ 

c 

They also give an estimate of convergence rate 

s < C2]6 - pOC 
— C\0. + C2p' 

This available analysis gives important guidance on designing preconditioners for MIN- 
RES method. That is, as long as we find an SPD operator M that satisfies the assumption 

(4.6) with constants c\ and c 2 independent of the physical and discretize parameters, we 
can use it to precondition MINRES and the convergence rate is uniform with respect 
to the discretization and physical parameters. According to [19], condition (4.6) means 
that operator H and AD 1 are norm-equivalent, and the corresponding preconditioners are 
known as norm-equivalent preconditioners. Next, we discuss how to apply such an idea to 
design practical and robust preconditioners for the MHD system. 

4.1.1. Application to the MHD System. From the discussion in the previous section, we 
know that based on the well-posedness of the problem, robust and effective norm-equivalent 
preconditioners can be obtained as long as the condition (4.6) is satisfied. In this section, 
we discuss how to use this framework to design preconditioners, which satisfy this con¬ 
dition, for the structure-preserving schemes (3.4) and (3.5). 

Preconditioner for auxiliary scheme (3.5). Following the proof of the well-posedness 
of the symmetric Picard linearization (3.2), we first design preconditioner for the auxiliary 
scheme with the stabilization term V ■ B, V ■ C). The operator form of (3.5) is 

( 4 . 7 ) 


( A\ 

-div* 

0 

^* \ 


( u \ 


( h \\ 

—div 

0 

0 

0 


V 


~g 

0 

0 

7 1 l3 + div*;i, 7 1 div) 

—a/iy 1 curl 


B 


-h 2 

V T 

0 

—acurl* fly 1 

SCTy J 


\e) 


\ h 3 / 


where 

A\u = H 1 u — Re~ 1 Au + k _ 1 div*divtt — scr r (u x B~ ) x B~ , Vu £ Vj u 
Tu = sa r u x B~, Vw £ V h , 

and hi = f, h 2 = ock~ l pj 1 B° , h 3 = 0, and g = 0. Note that, we have divh 2 = 0 and A\ = Hi- 
Based on the well-posedness of the scheme (see Theorem 3.1), A is an isomorphism 
from X h x Qi, to (X h x Q/,)*. Therefore, one simple choice of the preconditioner is the 
operator induced by the inner product (■, ■)-£ which satisfies (4.6) with constants c\ = c 2 - 
1, i.e. V = diag . This is because for the MHD system, || • ||^ is defined 

by (3.6). More precisely, operator V is a natural isomorphism from (X h x Q h )* to X h x Q h 
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with operator form 


(4.8) 


Mi 

o 

0 

Vo 


0 

kX 2 

0 

0 


0 0 \ 

0 0 

Ci{k~^ + div*fMdiv) 0 

0 sa r X 4 + fcacurl* ^p 1 curl / 


Based on previous analysis, we can directly conclude that V is a norm-equivalent pre¬ 
conditioner for A. 

A direct application of preconditioner V in practice can be expensive and time-consuming 
because we need to invert each diagonal block of V. In order to make the preconditioner 
more practical, we replace the diagonal blocks of V by their spectral equivalent SPD ap¬ 
proximations, i.e. 


(4.9) 


M = diag (Qi, Q 2/ Q 3 , Qi ), 


where 

(4.10) c 2 ,i ( QjX,x ) < [XLj 1 x / x S j < c 1/( - ( Qjx,x ), i = 1,2,3,4. 

In the implementation, we define Q\ by several MG cycles, Q 2 by simple iterative meth¬ 
ods such as Jacobi methods, and Q 3 , Q 4 by the well-known HX-preconditioner [15]. Then 
M satisfies the condition (4.6) with c\ = minjcpj, c\\, and c 2 = max{c^[, c^\, Cj 3 ,G 4 } 

Therefore, M is a norm-equivalent preconditioner for A. 

Preconditioner for scheme (3.4). Now we consider the original structure-preserving 
discretization without the stabilization term (/M V ■ B, V ■ C). Similarly, the operator form 
of (3.4) is 


(4.11) 


/ At 

—div* 

0 

T* \ 


h) 


( hl \ 

—div 

0 

0 

0 


P 


~g 

0 

0 

-afc _1 ^y 1 z 3 

— aj,i~ 1 cuil 


B 


-h 2 

V T 

0 

—acurl *fM 

S(7 r J 4 / 


\e) 


\ / 


Since V is a uniform preconditioner for A which has the stabilization term, it turns out 
that we can obtain a uniform preconditioner V for A by removing the stabilization term 
in V, i.e. 


(4.12) 


T> = diag 


nj — 1 nj — 1 
/L-i / T La / 


ock 1 }i r 1 X 3 


n -1 



Actually, using the fact that div curl = 0, we have 


a.(k 1 }i r M + div*/i r 1 div) (— oc}i r 1 curl) = —fccurl = ak 1 ji r 1 X 3 1 curl). 


Therefore, VA = VA and we conclude that V is a norm-equivalent preconditioner for 
A. Again, using V in practice can be expensive and time-consuming. We can replace the 
diagonal blocks of V by their spectral equivalent SPD approximations that satisfy (4.10), 
namely. 


(4.13) 


M = diag (Qi, Q 2/ Q 3 , Q 4 ). 


It is easy to see that At is a norm-equivalent preconditioner for A. 

A divergence-free preserving iterative process. One important feature of the structure 
preserving discretization is that it preserves the Gauss's law of magnetic field exactly on 
the discrete level. Therefore, when designing preconditioners and iterative methods, we 
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would like to inherit this property as well. We now discuss how to preserve this property 
mainly for scheme (3.4) and briefly for scheme (3.5) briefly. 

First, we consider a simple linear iteration based on operator V: 


( ul+1 ) 


(A 


( ei> 


/eA 



( h A 


(A 


p l+1 


p 1 

+ 

e 2 


e 2 

= v 


s 

- A 

p 1 


B ,+1 


B l 


e 3 


e 3 



-h 2 


B l 


\E m j 


\E l j 


\ e 4/ 


\ e 4/ 



\ h 3 y/ 


\E l ) 



It is easy to check that, if divB 1 = 0 and divh 2 = 0, then dive 3 = 0. Therefore, div£? ,+1 = 0. 
Namely, the linear iterative scheme based on B preserves the divergence-free condition 
at each iteration. 

Now we consider the preconditioned MINRES method with V as a preconditioner. The 
following theorem states that if the initial guess satisfies the divergence-free condition 
exactly, the solutions of all the iterations also satisfy this condition exactly. 

Theorem 4.1. If the initial guess x° = ( u° , p°, B°, E°) T , satisfies the divergence-free condition 
exactly, i.e. di vB° = 0, and d\wh 2 = 0, then all the iterations x l = (u 1 , p l , B 1 , E , ) T of the precon¬ 
ditioned MINRES method with preconditioner B satisfy the divergence-free condition exactly, i.e. 
di vB 1 = 0. 

Proof. According to the definition of preconditioned MINRES method with precondi¬ 
tioner V, we know that 

(4.15) x 1 ex Q + K\VA, r Q ), 
where 

tC l {VA, r°) = span{r°, VAr°, ( VA) 2 r°, ■ ■ ■, (: VA) l ~ l r °}, 

with r° - F — T>Ax°, r° = (r°, r\, r®, r°) T . Note that divr^ = 0 by its definition. 

Now, we denote (DA)' n r° by v m = (vf, vf, v™, v™) T for m = 0,1,2, ■ • ■, l — 1. Since v m = 
VAv'n- 1 and 

v 3 = (ftfc _1 ^,7 1 X3) _1 Oik^ p~p l l 3 V™~ 1 — ocpp 1 c urlv'f- 1 ^ 

(4.16) = -v" l ~ l - kcur\vf-\ 

divv’f - 0 if d'wvf- 1 = 0. Noticing that divv® = divr-g = 0, by mathematical induction, we 
have dwv'f = 0. 

Finally, due to the fact (4.15), x 1 is a linear combination of v m , m = 0,1 ,---,/ — 1 and then 
B 1 is a linear combination of vf, m = 0,1, • • •, l — 1. Because divv™ = 0, m = 0,1, • ■ ■, l — 1, 
we can conclude that divT? , = 0. □ 


Although using V as a preconditioner preserves the divergence-free condition exactly, 
in general, it is not true when we use operator M defined by (4.13) as preconditioner. In 

order to preserve this property, we need to use Q 3 = 

M, namely. 


oik p r I 3 in the preconditioner 


(4.17) 


M = diag ( Q 1 , Q 2 , ock 1 p r 1 1 3 


,0.4 


Remark 4.1. In terms of implementation, note that, for the magnetic field in v m = MAv m ~ x , we 
still have (4.16), which can he used to update v™ at each preconditioning step without inverting 
the mass matrix. 
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The following theorem states that if the initial guess satisfies the divergence-free con¬ 
dition exactly, then the solutions of all the iterations satisfy this condition exactly when 
M defined in (4.17) is used as a preconditioner. 

Theorem 4.2. If the initial guess x° = (u°, p°, B°, E°) T , satisfies the divergence-free condition 
exactly, i.e. di vB° = 0 and div /12 = 0, then all the iterations x l - (u ! , p 1 , B l , E 1 ) 7 of the pre¬ 
conditioned MINRES method with preconditioner M defined in (4.17) satisfy the divergence-free 
condition exactly, i.e. di vB l = 0. 

Proof. The proof is the same as the proof of Theorem 4.1 with V replaced by A4. □ 

Next we consider the linear system Ax = F . If we use V as the preconditioner, the 
solutions of all the iterations satisfy the divergence-free condition exactly if the initial 
guess is divergence-free. The argument is exactly the same with that of Theorem 4.1 
except that (4.16) is replaced by 

(4.18) vf = Uf 1 - ocp; 1 curb;”'" 1 ) = -vf- 1 - fccurk^" 1 . 

And if we use M, with Qj = Elf 1 = ^a(fc _ 1 //,T 1 + d'w' fif 1 d\v) j , as the preconditioner, 
i.e. 

(4.19) Af=diag(s 1 ,S 2 ,^ 3 - 1 ,S 4 ), 

all the iterations satisfy the divergence-free condition exactly as long as the initial guess 
is divergence-free. 

4.2. FOV-equivalent preconditioner. In this section, we recall the abstract framework 
for designing the FOV-equivalent preconditioners, following [19]. FOV-equivalent precon¬ 
ditioners are not necessary to be SPD, which makes it more general than norm-equivalent 
preconditioners. For example, we can design block triangular preconditioners, and use 
them to precondition the GMRES method. 

Consider the model problem (4.4), now A is not necessary to be symmetric. We use a 
general operator Me ■ H* —» H to denote the preconditioner. Based on the inner product 
(■, -) M ~i and the norm || • \\ M -i, we can estimate the convergence rate of the preconditioned 
GMRES. It is proved [11, 28] that if x m is the ///-iteration of GMRES method and x is the 
exact solution, then 

\\M c A(,x-x m )\\ M -' < _ f_\ m/2 

\\McA(x —x°)\\ M -\ ~ V V) ' 

where 

(420) < (x,MjcAx) m -i , \\M c Ax\\ m -i < ^ 

(x / x)^~ i II^IIm-i 

According to the theory, we conclude that as long as we find an operator Me arid a 
proper inner product (■, -) M i such that condition (4.20) is satisfied with constants 7 and 
F independent of the physical and discretize parameters. Me is a uniform preconditioner 
for the GMRES method. Such preconditioners is usually referred to as FOV-equivalent 
preconditioners. 
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To give a flavor of the analysis of FOV-equivalent preconditioners, we first demon¬ 
strate the analysis with a 2-by-2 block system. Application to the MHD system is dis¬ 
cussed later. Assume that Ax = F is in the form 


(4.21) 



-B* 

0 




where A\ is a SPD operator. Based on the partition of the system, we assume that a 
splitting of the Hilbert space H is Hi x H 2 such that x\ £ Hi, x 2 £ H 2 . Assume that 
the problem (4.21) is well-posed with respect to norm which is induced by M - 

diag ('H\,'H 2 )*~ 1 . And we further assume that = A\. Therefore, the well-posedness 
implies that there exists a constant £ > 0 , independent of physical and discretization 
parameters (depending on the problem) such that 


(4.22) 


inf sup 
x 2 £H 2 a-jgffj \\X\ 


(Bxi,x 2 ) 


* 2 ||h 2 


>£>o. 


Theorem 4.3. If the condition (4.22) holds, there exist constants 7 and T such that for all x f 0 , 
the operator A defined in (4.21) and the operator 


Me 




-1 


satisfy condition (4.20) with the norm || ■ \\ M -i induced by M = diag ( A\, Hi) 1 


Proof. By simple computation, we get 


1 Y2* 


McA = 


_ (Xi -AfyB 


0 HfyBAfyB* 


Then for any x = (x 2 , x 2 ) T , we have 

(x, M c Ax) m 1 = {xi - AfyB*x 2 , xi) Al + (BAf 1 B*x 2 , xi) 
= \\xi\\ 2 Ai -(B*x 2 ,xi)+\\B*x 2 \\ 2 A - 1 
> 11*1 - ll 33 ! IUj II ^** 2 11^-1 +1| B*a;2 ||^_i 


i; U 




where £1 = 11 m, £ 2 = || r/; 2 1| 1 . Since the matrix in the middle is SPD, there exists 


70 > 0 such that 


Moreover, 


(x,M c Ax) m -i > 7 o ( \\xi\\ z A +\\B*x 2 \\ 2 . 1 ) . 


1 12* 11 {Bx 1, xi) .. .. 

I B * 2 |U- 1 = sup > Q\x 2 \\ n2 , 

1 asie'Hi 11*1 Mi 


we get 

(x,M c Ax) m - 1 > 7oll*i|l^ 1 +7o? 2 ||*2||^ 2 > m in |7o, 7o£ 2 } (*, *)^-i. 


which leads to the lower bound 7 . The upper bound T follows directly from the bound¬ 
edness of each term. □ 
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As we mentioned before, applying Me defined in Theorem 4.3 as preconditioner can 
be expensive and time-consuming. Therefore, we replace the diagonal blocks by their 
spectral equivalent SPD approximations. The following theorem states that under certain 
assumptions, such preconditioner is still robust. 


Theorem 4.4. If the condition (4.22) holds, there exist constants 7 and T such that for all x f 0 , 
the operator A defined in (4.21) and the operator 


Me 

satisfy condition (4.20) with the norm || 


_ 21 


-i 



B Q~fy) 

induced by M - diag (Q lf Q 2 ) provided that 


(1) C 2 ,i ( QjX,x ) < (fH i 1 x,x s j < c u ( QiX,x ), i = 1 or 2, 

(2) ||Ii - QiA\\\^< p, with 0 < p < 1 . 


Proof By simple computation, we get 

MeA-( M ~ Q ^\. 

- Q 1 A 1 ) Q 2 BQ 1 B*J 

Then for any x = {x\, x 2 ) J , we have 

(x,M c Ax) m -! = \\x 1 \\ 2 Ai -(B*x 2 ,x 1 ) + {B(I 1 - Q 1 A 1 )x 1 ,x 1 )+ \\B*x 2 \\ 2 Qi 
= \\xi\\ z Ai -(QiA 1 x 1 ,B*x 2 )+ \\B*x 2 \\ 2 Qi . 

As ||Xi - Q\A\ < p implies that 


(l-p)(*l,*l)^-i < (Xl,X 1 ) Ql < (1 + p)(x 1 ,x 1 ) A 1, 

(1+ l o) _ 1 (a; 1 ,£c 1 )^ 1 <{x lr Xi) Q -i < (1 - pT l {x lr xf) Al , 

we have 


— {Q\A\X\, B*X 2 ) < Ml®l||Q 1 ||iS*®2||Q 1 < (l+PlMl^lll^-l ||^** 2 || fii 
= ( 1 +P)\\xi\\a 1 \\B*x 2 \\q 1 , 


Therefore, 


(x,M c Ax) m -i > \\x 1 \\ 2 A ^-(l+p)\\x 1 \\ Al \\B*x 2 \\ Ql + \\B*x 2 \\ 2 Qi 

(SlY ( 1 -Q+p)/2\ (&\ 

\S 2 ) \-V + P)/2 1 ) W' 


where = \\xi 1^, £ 2 = \\B''x 2 \\q v We can verify that the matrix in the middle is SPD 
when 0 < p < 1 . Therefore, there exists a constant 70 > 0 such that 

{x,M c Ax) m -i > 70 [\\xi\\\+\\B*x 2 \\ 2 Ql ^ > 7 o(l - P)\\xi || g-i+ 7 o(l - pK 2 \\x2\\h 2 

> min|7o(l -p), 7 o(l ~ pK 2 cf}} (x,x) M -i, 

which leads to the lower bound 7 . The upper bound T follows directly from the fact that 
each term is bounded. □ 


We comment that the second assumption in Theorem 4.4 is reasonable as in practice 
we can achieve it by performing one or several steps of V-cycle multigrid method. 
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4.2.1. Application to MHD system. In this section, we discuss how to design FOV-equivalent 
preconditioners for the structure-preserving discretization (3.4). Instead of giving details 
of the FOV-equivalent preconditioners for the auxiliary problem (3.5), we comment that 
similar preconditioners exist, and theoretical results follow by similar argument. 

The operator form of the mixed formulation (3.4) is 


(Ax 

—div* 

0 

E* \ 


/ M \ 


( hl \ 

div 

0 

0 

0 


V 


g 

0 

0 

«A: _1 pT 1 l3 

app'curl 


B 


h 2 

\r 

0 

—ftcurFpp 1 

scr r X 4 / 


\e) 


Us/ 


Here, we multiply -1 to the second and the third equations because we do not require 
the system to be symmetric any more and such a small modification makes our exposition 
more clear. Moreover, we still have T~L\- A\. 

Based on the well-posedness of scheme (3.4), it is natural to propose the following 
block lower triangular preconditioner Me- 


(4.23) 


M c 


(A x 0 0 0 \ 

div kX 2 0 0 

0 0 0 
\ E 0 — acurFpy 1 XL 4 ) 


Note that the diagonal blocks of M c ] are the same as those of V 1 defined by (4.12). 
Thus, naturally we choose inner product (■, ■)% and norm ||- || w (3.6) to verify the condi¬ 
tion (4.20). We prove that operator Me is a robust FOV-equivalent preconditioner in the 
following theorem. 


Theorem 4.5. If k < ko, there exist constants 7 and T that are independent of the mesh size h, 
time step size k, mid physical parameters Rm, s, p r , and <x r , such that for all x f 0, the condition 
(4.20) holds with (-, -) M -x = (•, -) H . 


Proof. By simple computation, we get 


M C A = 


fa 

—A 1 1 div* 

0 

A^E* 

0 

A: _1 divHp 1 div* 

0 

-k^dwA^E 

0 

0 

Ts 

/ecu r 1 

\o 

H^EA^d iv* 

0 

h a 1 s 4 


where 

S 4 = sa r X 4 + cdccurFp^curl — EAf^E*. 

Then for any x = (u , p, B, E) T , we have 

(; x,McAx)'h = (u,u) j4i — (u, div*p) + (u, E* E) + (p, div7lp 1 div*p) 

- (p, di vA^E*E) + (B, B) n 3 + cc(B,cur\E) jh l 
+ (E, EAf l dw*p) + (E, S 4 E) 

> ll«ll?t 1 -|l M IU 1 l|div*p||^ r i-||tx|U 1 ||J c '*i;||_ Ar i + ||div*p||^_ 1 + ||B||| 1 

- \\B\\h aV^Ucurl^H i+s\\E\\l+ock\\cur\E\\ 2 4 -\\E*E\\ 2 

rr p r 
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Since 


,*-**,„ (E*E,v) (E,Tv) ^ ^\\E\\ ar ^\\v x B- 

\J- E\\a-i = sup -= sup - < sup —-— 7 i—^- 

1 v=/0 MUi v=/0 MUi vJO 


\ v \\A 1 


< ~^V~s\\E\\ ar , 


we get 


Vs, 


{x,M c Ax) n > ||w||ji 1 -||w|U 1 ||div*p|| iA -i--^||w|U 1 ||£|| (ri .+||div*p||^_ 1 +||B ||^ 3 
- \\B\\ H3 Voik\\cuv\E\\ 1 + S -\\E\\l+ak\\cuilE\\ 21 t 




T 

/ 1 -1/2 

0 -1/2V2 

0 > 


52 


-1/2 1 

0 0 

0 

= 

53 


0 0 

1 0 

-1/2 


54 


- 1 / 2 V 2 0 

0 1/2 

0 




v 0 0 

-1/2 0 

1 

u \\au 5i = 

\div*p\\ A -i, 53 = \\E\\n 3 , 54 = a/sII-EIU 

and £5 


(ti\ 

52 

53 

54 

\5sJ 


'Mr 


is easy to verify that the matrix in the middle is SPD. Therefore, there exists a constant 
7o > 0 such that 

(x,M c Ax) n >70 ^||u||^ li + ||div*p||^_ 1 + ||S ||^ 3 + ||£;||^ 4 ^ > min [ 70 , 70 ^} (x,x) M -i. 

The last inequality comes from the fact that 

11 * n (div*p, v) (p, divu) 

||div p\\ A 1= sup = sup > 5 Mh 2 - 

1 170 \\ v \\Ai 170 II^IUi 

The above estimate leads to the lower bound 7 . On the other hand, the upper bound T 
can be obtained directly by the boundedness of each term. □ 


To reduce the time and computation cost of Me, we replace its diagonal blocks by 
their spectral equivalent SPD approximations except that of B. The implementation is 
the same as that in remark 4.1. Such modification gives rise to the following operator 


(4.24) 


Me 


/ 2r 1 

0 

0 

° ^ 

div 

S 2 1 

0 

0 

0 

0 

X 3 

0 

V A 

0 

—acurP^y 1 

S 4 " 1 / 


We note that, for the magnetic field, we can apply the same approach with (4.16) at each 
preconditioning step in order to preserve divergent-free condition exactly on the discrete 
level. And we can prove a similar conclusion to that in Theorem 4.5 for Me using a 
different inner product in the analysis. 


Theorem 4.6. If k < ko and the condition (4.10) holds, then there exist constants 7 and Y that 
are independent of the mesh size h, tune step size k, and physical parameters Rm, s, y r , and a r , 
such that for all x f 0 , the condition (4.20) holds with the inner product (x,y) M -i induced by 
M = diag (^Qi, 0.2,'Hf 1 , Q 4 ) provided that \\Xi - Q 1 A 1 II 7 < p, with 0 < p < 0.289. 
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Proof. By calculation, it can be seen that 


M C A = 


( QiAr 
Q2div(Xj — Q\A\) 
0 

V i - Q.At) 


-Qtdiv* 

Q.2&2 

0 

Q^Qidiv* 


0 

0 

X 3 

0 


QiJ r * \ 

-Q 2 divQ 1 T* 
/ecu r I 

QiSi J 


where S 2 = divQidiv*, <S 4 = so>X 4 + a/cciirl '/y 'curl - EQ\E*. Therefore, for any x = 
( u , p, B, E) T , 


{x, McAx) m -i = (u, A\u) - (it, div*p) + (u , F*E) + (p, div(X 4 - Q-lA^u) + (p, S 2 p ) 

- (p,divQ 1 X'*.E) + (B,B) n3 + (B,kcur\E) H3 + (E, T(fL\ - SiTliM 
+ (JS, J^Qidiv’-p) + (J5, £ 4 #) 

> llwll^-llwlUill-^-SlUj-i-MiulleilldivVIlej + lldivVIl^ 

+ II-®IIh 3 —ll-^IIWaV^llcurl-Ell^-i—p||u.||_4 1 1| J r *£ ! ||^-i 
+ S || J E||2 r+a fc||curl£;||J rl -||^£;|| 2 Qi . 

Since ||Xi - Qi-4i|Uj< p implies 


(1 -p)(v,v) A -i < (v,v) Ql < (1 + p)(v,v) A -i, 

(1 + p)~\v,v) Al < {v,v) Q -i < (1 -prHv,v) Al , 


and \\E*E\\ a -i< 


1 

71 


Vs\\E\U 


(it is shown in the proof of Theorem 4.5), we get 


_ 

(x,M c Ax) m -i > H M ll3l 1 -7|ll“IU 1 ll^llcr r -(l + p)||w|U 1 ||div>|| Cl +||div>||| i 


+ \\B\\n-\\B\\ H3 VKk\\curlE\\ i-p^L\\u\\ Al \\E 


IW 3 

+ li y 0 s||£;|| 2 r +afc||curl£;|| 2 _ 1 


72' 


& 

£3 

& 

V&7 


/ 


—(1 + P )/2 
0 

-(l+p)/2x/2 

0 


—(1 + p)/2 0 ~(l + p)/2y/2 0 \ 

10 0 0 
0 1 0 - 1/2 

0 0 (1 — p )/2 0 

0 - 1/2 0 1 ) 


/SA 

£2 

£3 

£4 

V& / 


where £1 = IMI^, £> = ||div*p|| Ql , £3 = ||B|| W3 , £4 - Vs||^IU and £5 - VafcUcurlSl^-i. we 
can verify that the matrix in the middle is SPD when 0 < p < 0.289. Therefore, there exists 
a constant 70 > 0 such that 


(x,M c Ax) m -i > 70 (^||«|l3i 1 +(l-|0)l|div>||5 l _ 1 +||B||^ 3 +||X;||^ 4 ^ 

> to (llull^+a - p)£ 2 IIpIIm 2 +I|£|I«3+II-e|I«4) 

> min |7c(l - p), 7o(l - p)^ 2 c~^, 7oc/]} (x,x) M 4 , 

which leads to the lower bound 7. On the other hand, the upper bound T can be obtained 
directly from the boundedness of each term. □ 
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Remark 4.2. Condition \\li - Qi^ilUj^ p, with 0 < p < 0.289, means that we should solve the 
velocity block accurately enough. We comment that this condition can be relaxed by introducing 
parameters p 4 and p 2 in front of the diagonal block Q 2 and Q 4 , respectively, In that case, the proof 
is very similar to that of the Theorem 4.6 by choosing different numbers in the Young's inequality 
and picking appropriate p\ and pi- Thus, we omit the proof. 


5. Generalizations to other discretizations 


In this section, we generalize the robust preconditioners based on the well-posedness 
to other discretizations of incompressible MHD systems that have been developed in the 
literature. As examples, we will give detailed discussions on the discretization proposed 
by Gunzburger et al. [14] for the stationary incompressible MHD system and that for 
non-stationary incompressible MHD system given in [3, 24], Generalizations to other 
discretizations are similar. 


5.1. H 1 discretization for a stationary incompressible MHD system. The stationary in¬ 
compressible MHD equations discussed in [14] is: 


(5.1) 


1 

u ■ Vit — —Am + Vp — a(V x B) x B = f, 
V ■ u = 0, 

1 

— V x (it x B) + -V xVxB = 0, 

v ' Rm 

, V B = 0, 


with homogeneous boundary conditions 


u = 0, n ■ B = 0, n x (V x B) - 0. 


The finite element space for ( u, p) is still V h x Q h . However, due to the boundary condition 
of B, a special Sobolev space is used, 

Hl(Cl) 3 = {u 6 HW, n ■ u |an= 0 } , 

and norm for this space is 


(5.2) 


u 


2 0C 

n ’ 1= Rm 


V • m|| 2 +|| V x u 


v« e Hjfn,) 3 . 


Then the finite element space for B is Hj l h (n) 3 c HjfCl) 3 . 

The full discretization based on symmetric Picard linearization of (5.1) is: find (uj,, Bj,, pf) e 
V h x H 3 h (Cl) 3 x Q h such that for any (v h ,C h ,q h ) e V h x Hj l h (Ci) 3 x Q,„ 

{ d(uf,uf, v h ) + ^(Vui, V« h ) - a(V x B,„ Bf x v h ) - (p h , V ■ v h ) = (f h , v h ), 

-oc(u h x B~, V x C h ) + ^ [( v x B lu V x C h ) + (V ■ B h , V ■ C h )] = 0, 

(V ■ rt; ; , qh) — 0, 


where u h and B h stand for solutions from last nonlinear iteration step. The operator 
form of this system is (we flip the signs of equations to make the system symmetric): 


( -Re- 1 A 

—div 

y yaB- x curl) 


—div* 

0 

0 


ctB 


x curl 
0 


\ 


— —— (div*div + curl*curl) 
Rm J 



C\x — F\ =L 
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In [14], it is proved that the above discretization is well-posed with respect to the 
canonical norms of Vj lr Q/, and the norm defined by (5.2) for ff,' ;i (Q) 3 . This implies ro¬ 
bust norm-equivalent preconditioner as follows: 




(-Re' 1 A 0 0 

0 X 2 0 

0 0 —— (div*div + curPcurl) 

Rm v ; 


-l 


As mentioned before, we can design FOV-equivalent preconditioner for the original non- 
symmetric system, which is given by 


■M-c, i 


( -Re -1 A 0 0 \ 1 

div X 2 0 

(oiB~ X curl) 0 —— (div*div + curPcurl)^ 


Use the similar argument in previous sections, we can show that these preconditioner is 
uniform with respect to h. The first diagonal block is Poisson-like, for which multigrid 
methods are efficient. Jacobi iterative method or other simple iterative methods are effi¬ 
cient for the second block. And the third diagonal block is a vector Laplacian, for which 
multigrid methods are effective as well. 


5.2. H(curl) discretization for a non-stationary incompressible MHD system. The math¬ 
ematical model discussed in [3, 24] is 


(5.3) 


du 1 

— + u ■ Vm — —Ait + Vp — a(V x B) x B — f , 
V ■ u = 0, 

SB 1 

— -V x (u x B) + — —V x V x B — Vr = 0, 


dt 

V B = 0, 


Rm 


with boundary conditions 


u = 0, n x B = 0. 

The full discretization based on symmetric Picard linearization of (5.3) is: find (u\\, p n h , B r f [, r)]) e 
Vh x Qh x V£ x Hi such that for any (v h , q h , C h , l h ) e V h x Q h x Vfi x H\, 

k- 1 (< - m"- 1 , Vhj + d(u\ l -\ u;-\v h ) + ^ (Vmj, Vv h ) 

-a ((V x B h ) x Bl~\v h ) - (pi V • v h ) = (fZ,v h ), 
ak- 1 (b% - B n h -\ C h ) - a ( u n h X B\\~ l , V x C H ) + ^(V x B\\, V x C, ; ) 

~oc{C h , Vr h ) = 0, 

(V ■ Uj ir q h ) = 0, 

*{Bl,\7l h ) = 0. 


(5.4) 
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The operator form of discretization (5.4) is (we flip the signs of equations to make the 
system symmetric): 

(k-% - Re- 1 A — div* 


C 2 x = F 


— div 

(kB- x curl)* 
0 


olB x curl 
0 

— 3 + —— curl* curl ^ agrad 


(agrad)* 


0 


\ 

/u\ 


( h l \ 


V 


gi 


B 


—h 2 

/ 

\r) 


\ g2 J 


Based on the analysis in the literature [3, 24], we can prove that discretization (5.4) is 
well-posed with respect to weighted norms in the finite element spaces Vj u Qh, Vfi, and 
Hj t . Therefore, robust norm-equivalent preconditioner for this system is: 

-l 


B 2 = 


(k 1 X\—Re 1 A + rdiv*div 0 
0 kX 2 

0 

0 0 


0 0 

0 0 

0 3 + —— curl* curl 0 

d Rm 


\ 


0 ka~ 1 (I A -A)J 

And an FOV-equivalent preconditioner for the original non-symmetric system is: 


■Me, 2 = 


(k 1 T\—Re 1 A + rdiv*div 0 
div kX 2 


—(aB X curl)* 0 ctk 1 X 2 + curl*curl 


\ 


0 


0 


(agrad)* 


k<x. x (X 4 — A)J 


When k is sufficiently small, we can show these preconditioner is robust with respect 
to discretization parameters. Multigrid methods are effective to precondition the first 
diagonal block. Jacobi method or other simple iterative methods are efficient for the 
second diagonal block. HX preconditioner is powerful for the third diagonal block. And 
the fourth diagonal block corresponds to a Poisson problem for which multigrid methods 
are powerful as well. 


6. Numerical experiments 

In this section, we carry out numerical experiments for both 2D and 3D incompress¬ 
ible MHD systems to demonstrate the robustness of the preconditioners introduced in 
previous sections. As mentioned before, we mainly focus on the structure-preserving 
discretization without the stabilization term, namely, scheme (3.4). 

In both 2D and 3D implementation, we use P 2 -Po to discretize the velocity and the 
pressure, and use the lowest order elements to discretize the electric and the magnetic 
field. We implement the structure-preserving discretization based on the FEniCs [18] and 
the block preconditioners based on the FASP package [34], We run all the experiments 
on a Dell workstation with 12 GB total memory. 

6.1. Implementation of block preconditioner. In this section, we discuss implementa¬ 
tion details of the proposed block preconditioners when they are used to precondition 
the Krylov methods such as MINRES, GMRES, and flexible GMRES (FGMRES). As we 
know, inverting block diagonal or triangular preconditioners ends up with inverting di¬ 
agonal blocks. Therefore, the main implementation issue is how to invert those diagonal 
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blocks. When inverting each diagonal block exactly, we call direct solvers implemented 
in the UMFPack package [10]. While when inverting each diagonal block approximately, 
we mainly call preconditioned conjugate gradient (PCG) method with a relative large tol¬ 
erance, since each block in our preconditioners is SPD. For example, the tolerance of PCG 
in terms of I 2 norm of the relative residual is 10~ 3 , which is sufficient to achieve (4.10). 

It is well-known that effective preconditioners for PCG method are necessary in order 
to achieve overall robustness of the preconditioners and reduce the computational cost 
as well. In both 2D and 3D, since the velocity block is a Poisson-like problem, multigrid 
methods are good candidates. In our experiments, we use AMG method implemented 
in FASP package. For the pressure, the corresponding block is a mass matrix of finite 
element of Q;„ we use Jacobi preconditioner for it. And for the magnetic block, we up¬ 
date it using (4.16) without inverting the diagonal block explicitly. This not only ensures 
the exactness of the divergence-free condition for the magnetic field but also saves the 
computation time. 

However, the implementations of the diagonal block corresponding to the electric field 
are different in 2D and 3D. This is because the vector electric field degenerates to a scalar 
in 2D. In 3D, the diagonal block for the electric field E corresponds to the following PDE 

a 2 curl curlD + sE = f. 


Therefore, HX preconditioner [15] implemented in FASP package is a powerful precon¬ 
ditioner for it. However, in 2D, the diagonal block of the electric field £ corresponds 
to « 2 rot curlE + sE - f on the continuous level. As rot curl/; = — An, this PDE is actually 
—ocAE + sE = f, which is a Poisson-like problem. Obviously, robust preconditioner for the 
above problem depends on the ratio between the parameter a and s. When a is relatively 
large, — A dominates and we use AMG as the preconditioner. When a is relatively small, 
the lower order term dominates and we use Jacobi as the preconditioner. 


6.2. Convergence tests. We first perform convergence tests to verify the correctness of 
our code. For the sake of accuracy, we use 2-step BDF for the temporal discretization. 
All the physical parameters are set to be 1 in the convergence tests. The initial condi¬ 
tions, boundary conditions, and right hand sides are computed based on the given exact 
solutions. The tolerance for the preconditioned Krylov methods is 10” 10 . 


6.2.1. 2D Convergence Test. The exact solutions chosen for the 2D convergence test are: 


u=r cos! 'j, bJ 0 '' 

\ 0 \ sin f cos x , 


—xcosy, £ = sinx. 
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(a) Error versus mesh size h ( (b) Error versus time step size 

k = 0.01 and t = 0.1) k (h = 1/32 and t = 0.8) 


(c) ||divB|| (A: = 0.01) 


Figure 6.1. Numerical results of 2D convergence test. 


Based on the results shown in Figure 6.1 (a) and (b), we can see that both spatial and 
temporal errors converge in the optimal order, i.e., 0(h) and 0(k 2 ). We also plot |div£?| 
in Figure 6.1 (c) which is about 10 12 and confirms the exactness of the divergence-free 
condition due to the structure-preserving discretization. One comment is that ||div.B|| 
increases as the time evolves. This is due to the accumulation of the round-off error 
during the computation. We can reduce such accumulation simply by updating Bjf with 

(6.1) B'fl =Bj-Vx 

at the end of each time step. (6.1) is a mathematical equivalent formula to the original 
system. Therefore, we can preserve the divergence-free condition exactly on the discrete 
level. 


6.2.2. 3D Convergence Test. The exact solution chosen for this test is: 


u 




p = — xcos y. 




(a) Error versus mesh size h (b) Error versus time step size 
(k = 0.01 and f = 0.1) k (h = 1/12 and t = 1) 


(c) ||divB|| (A: = 0.01) 


Figure 6.2. Numerical results of 3D convergence test. 


Based on the results shown in Figure 6.2, we can see that both spatial and temporal 
errors converge in the optimal order and the divergence-free condition is preserved. 
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6.3. Benchmark problem: lid-driven cavity. Next we consider the lid-driven cavity prob¬ 
lem, which is a well-known benchmark problem for CFD simulation [29]. A constant 
magnetic field Bq is applied and coupled with fluids. When the external imposed mag¬ 
netic field is zero, the problem becomes the classical hydrodynamic lid-driven cavity 
problem. Here, we assume that the cavity is a unit square in 2D and a unit cubic in 3D. 


6.4. 2D lid-driven cavity. For the 2D case, the fluid movement in the cavity is induced by 
the imposed boundary condition u = u 0 on the boundary y = 1 (figure 6.3). The constant 
background magnetic field is B t] = (0,1) T . The velocity field perturbs to the magnetic 
fields, and the Lorentz force acts on the fluid and changes the motion of the fluid flow. 


x 


Figure 6.3. Geometry of lid driven cavity 

Again, we use symmetric Picard iteration in this benchmark test. Parameters for this 
test are: 

• Time step size k = 0.01, time interval is [0,10.0]. 

• Uniform grid is used, with h = 1/100. Number of elements is 2 x 10 4 . 

• The tolerance of preconditioned Krylov methods is 10~ 8 and the tolerance of non¬ 
linear iteration is 10 6 . 

• }i r = 1, cr r = 1, and s = 1. 

We perform numerical tests of different Reynolds number (Re) and magnetic Reynolds 
number ( Rm). In Figure 6.4, we plot stream lines of velocity fields in different cases. We 
can see the vortex appears as the Reynolds number increases. 



400 1 400 


FIGURE 6.4. Stream lines of velocity field. 


In Figure 6.5, we plot the distribution of the total magnetic fields for different Re and 
Rm. When Rm is relatively small, the total magnetic field is almost the same as the con¬ 
stant background field, while when Rm is relatively large, the distribution of total mag¬ 
netic field is attached to the motion of the fluid. 
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1 400 


FIGURE 6.5. Distribution of total magnetic field. 


Next we demonstrate the performance of the proposed preconditioners V (4.12) and 
Me (4.23) in Table 6.1. Note that, these two preconditioners requires to solve each diag¬ 
onal block exactly. The numbers in the table are the number of iterations for precondi¬ 
tioned Krylov methods. 



V (MINRES) 

V (GMRES) 

M c (GMRES) 

Re = 1, Rm = 1 

32 

27 

7 

Re = 1, Rm = 400 

33 

26 

6 

Re = 400, Rm = 1 

25 

21 

5 

Re = 400, Rm = 400 

27 

23 

5 


TABLE 6.1. Performance of solver (2D): s = 1 


Besides the performance of solver, we also concern about the value of k 0 , defined by 
(3.7). So we compute the values of ko for each time step and nonlinear step, and list the 
minimum value of ko of all the iterations in table 6.2. Based on the values of ko, we can 
conclude that the refinement of the mesh does not influence the restriction on time step 
size significantly. 



Re = 1, Rm = 1 

Re = 1, Rm = 400 

Re = 400, Rm = 1 

Re = 400, Rm = 400 

h = 50 

0.12 

0.00084 

0.12 

0.0017 

tr 

ii 

o 

o 

0.12 

0.00065 

0.12 

0.0012 

h = 200 

0.12 

0.00057 

0.12 

0.00098 


Table 6.2. Restriction on time step size: values of ko 


In order to further investigate the robustness of the proposed preconditioners with 
respect to the time step size k and the mesh size h, we list the number of iterations of 
the Krylov methods with different block preconditioners in Table 6.3, 6.4, 6.5, 6.6 and 6.7. 
The tolerance of the preconditioned Krylov methods is 10 6 . When the diagonal blocks 
are solved approximately by the PCG method, the preconditioners are actually changing 
in our implementation. Therefore, we use FGMRES method in this case to improve the 
robustness. Based on the results shown in the tables, we can conclude that our block 
preconditioners are effective and robust. 
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(a) Re 

= 1 , 

Rm = 

1 

(b) Re 

= 1, 

Rm = 

400 

(c) Re 

O 

O 

II 

Rm 

= 1 

(d) Re = 

400, 

Rm = 

400 

/z 

Af^\ 

1/32 

1/64 

l/l28 

h 

A 

1/32 

l /(A 

!/l28 

h 

A 

1/32 

i/64 

1/128 

h 

A£^\ 

!/32 

!/64 

l/l28 

0.02 

23 

23 

22 

0.02 

19 

23 

23 

0.02 

15 

15 

15 

0.02 

14 

16 

17 

0.01 

24 

24 

22 

0.01 

17 

19 

21 

0.01 

13 

15 

15 

0.01 

10 

14 

16 

0.005 

24 

23 

22 

0.005 

16 

17 

19 

0.005 

13 

13 

15 

0.005 

10 

12 

15 

0.0025 

23 

23 

22 

0.0025 

14 

16 

17 

0.0025 

13 

13 

13 

0.0025 

9 

10 

13 


Table 6.3. Block diagonal preconditioner V (4.12) for MINRES method 
(diagonal blocks are solved exactly) 


(a) Re 

= 1 , 

Rm = 

1 

(b) Re 

= 1 , 

R»z = 

400 

(c) Re 

O 

O 

II 

Rm 

= 1 

(d) Re = 

400, 

Rm = 

400 

/z 

Af^\ 

1/32 

1/64 

l/l28 

h 

1/32 

1/64 

l/l28 

/z 

A£^\ 

1/32 

i/64 

1/128 

A£^\ 

1/32 

1/64 

l/l28 

0.02 

20 

18 

17 

0.02 

16 

17 

15 

0.02 

13 

11 

11 

0.02 

8 

11 

9 

0.01 

19 

18 

18 

0.01 

14 

15 

15 

0.01 

11 

13 

12 

0.01 

7 

7 

7 

0.005 

21 

20 

18 

0.005 

14 

14 

15 

0.005 

9 

12 

12 

0.005 

6 

7 

7 

0.0025 

20 

20 

18 

0.0025 

12 

14 

13 

0.0025 

10 

9 

10 

0.0025 

8 

7 

7 


Table 6.4. Block diagonal preconditioner V (4.12) for FGMRES method 
(diagonal blocks are solved exactly) 


(a) Re 

= 1, 

Rm = 

1 

(b) Re 

= 1, 

Rm = 

400 

(c) Re 

O 

O 

II 

Rm 

= 1 

(d) Re = 

400, 

Rm = 

400 


1/32 

1/64 

l/l28 

/l 

1/32 

1/64 

l/l28 


1/32 


1/128 

h 

1/32 

1/64 

l/l28 

A£^\ 


A£^\ 

1/64 


0.02 

6 

5 

5 

0.02 

5 

5 

5 

0.02 

4 

4 

4 

0.02 

3 

3 

3 

0.01 

5 

5 

5 

0.01 

5 

5 

5 

0.01 

3 

4 

4 

0.01 

3 

3 

3 

0.005 

5 

5 

5 

0.005 

5 

5 

5 

0.005 

3 

3 

3 

0.005 

3 

3 

3 

0.0025 

5 

5 

5 

0.0025 

5 

5 

5 

0.0025 

4 

3 

3 

0.0025 

4 

3 

3 


Table 6.5. Block lower triangular preconditioner Me (4.23) for FGMRES 
method (diagonal blocks are solved exactly) 


(a) Re 

= 1 , 

Rm = 

1 

(b) Re 

= 1, 

Rm = 

400 

/z 

A£^\ 

1/32 

1/64 

l/l28 

k 

A f 

1/32 

1/64 

l/l28 

0.02 

21 

25 

43 

0.02 

16 

22 

43 

0.01 

22 

29 

51 

0.01 

16 

24 

49 

0.005 

24 

34 

60 

0.005 

16 

25 

57 

0.0025 

24 

38 

70 

0.0025 

15 

27 

60 


(c) Re 

O 

O 

II 

Rm 

= 1 

(d) Re = 

400, 

Rm = 

400 

/z 

A£^^^ 

1/32 

1/64 

l/l28 

k 

At 

1/32 

1/64 

l/l28 

0.02 

16 

24 

46 

0.02 

13 

23 

46 

0.01 

16 

26 

47 

0.01 

12 

21 

47 

0.005 

18 

24 

49 

0.005 

12 

20 

47 

0.0025 

20 

25 

49 

0.0025 

14 

21 

47 


TABLE 6.6. Block diagonal preconditioner M (4.13) for FGMRES method 
(diagonal blocks are solved approximately) 


(a) Re 

= 1 , 

Rm = 

1 

(b) Re 

= 1, 

Rm = 

400 

(c) Re 

= 400, 

Rm 

= 1 

(d) Re = 

400, 

Rm = 

400 

/z 

Af^\ 

1/32 

1/64 

l/l28 

/l 

1/32 

1/64 

l/l28 

/z 

A£^^^ 

1/32 

1/64 

l/l28 

h 

A 

1/32 

1/64 

l/l28 

0.02 

6 

9 

15 

0.02 

6 

9 

15 

0.02 

6 

9 

16 

0.02 

6 

9 

16 

0.01 

6 

9 

17 

0.01 

6 

9 

17 

0.01 

6 

10 

17 

0.01 

6 

10 

17 

0.005 

6 

10 

20 

0.005 

6 

10 

20 

0.005 

6 

10 

17 

0.005 

6 

10 

17 

0.0025 

6 

11 

22 

0.0025 

6 

11 

22 

0.0025 

7 

11 

18 

0.0025 

7 

10 

18 


Table 6.7. Block lower triangular preconditioner Me (4.24) for FGMRES 
method (diagonal blocks are solved approximately) 
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6.5. 3D lid-driven cavity. In 3D case, the fluid movement in the cavity is induced by 
the imposed boundary condition u = uq - (1,0,0) r on the plane y - 1. The constant 
background magnetic field is B 0 = (0,1,0) T . 


y 


B 0 


u 0 




FIGURE 6.6. Geometry of lid driven cavity 

We also use symmetric Picard iteration to solve this 3D benchmark problem. Since the 
motion of the fluids is uniform along the z-axis, the sectional view along z-axis of the 
motion of the fluids and the distribution of the total magnetic field are similar to those of 
the 2D benchmark. Therefore, we only present the performance of our preconditioners 
in the 3D benchmark. Parameters for the 3D case are: 

• Time step size At = 0.01. Time interval [0,1.0]. 

• Uniform grid with mesh size h = 1/12. Number of elements is 10368. 

• The tolerance of the preconditioned Krylov methods is 10” 8 and the tolerance of 
nonlinear iteration is 10 ~ 6 . 

• }i T = 1, a r = 1, s = 1. 

For 3D, we use coarser mesh than the 2D tests (but compatible problem size) due to 
the limited memory capacity of our workstation. We list the number of iterations of the 
Krylov solvers with preconditioners V (4.12) and Me (4.23) in Table 6.8. Note that, we 
need to invert each diagonal block exactly for these preconditioners. 



V (FGMRES) 

M c (GMRES) 

Re = 1, Rm = 1 

160 

36 

Re = 1, Rm = 400 

86 

35 

Re = 400, Rm = 1 

54 

16 

Re = 400, Rm = 400 

34 

16 


TABLE 6.8. Performance of solver (3D): s = 1 


In order to further demonstrate the robustness of the preconditioner with respect to 
the time step size k and the mesh size h, we list the number of iterations of the Krylov 
methods with different block preconditioners in Table 6.9, 6.10, 6.11, 6.12, and 6.13. The 
tolerance of the preconditioned Krylov methods is 10 _6 . * in the tables means that our 
workstation is out of memory during the test. This usually happens when we call direct 
solvers. In 3D, diagonal blocks get bigger when the mesh size decreases, so we cannot 
call direct solvers due to the limitation of our workstation. When solving each diagonal 
block approximately, we set the tolerance of PCG 10 3 . Based on the test results, we can 
conclude that our block preconditioners are effective and robust. 
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(a) Re 

= 1, 

Rm = 

1 

(b) Re = 

1, Rm = 

400 

(c) Re = 

400, 

Rm 

= 1 

(d) Re = 400, 

Rm 

= 400 

h 

A 

V 4 
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1/16 

h 
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1/4 

1/8 

1/16 

h 

A 
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V» 

1/16 
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A 
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Vs 
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65 

79 

* 
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44 

56 

* 

0.02 

45 

29 

* 

0.02 

32 

21 

* 

0.01 

59 

72 

* 

0.01 

38 

47 

* 

0.01 

48 

32 

* 

0.01 

31 

23 

* 

0.005 

54 

67 

* 

0.005 

34 

44 

* 

0.005 

49 

36 

* 

0.005 

32 

23 

* 

0.0025 

42 

61 

* 

0.0025 

28 

40 

* 

0.0025 

49 

40 

* 

0.0025 

32 

25 
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Table 6.9. Block diagonal preconditioner V (4.12) for MINRES method 
(diagonal blocks are solved exactly) 


(a) Re 

= 1, 

Rm = 

i 

(b) Re = 

1, Rm = 

400 

(c) Re = 

400, 

Rm 

= 1 

(d) Re = 400, 

Rm 

= 400 

h 

At^\^ 

1/4 
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1/16 

h 

A 

1/4 

1/8 

!/ l 6 

h 
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'A 
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54 

62 

* 

0.02 

34 

40 

* 

0.02 

30 

19 

* 
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21 

13 

* 

0.01 

52 

49 

* 
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32 

28 

* 

0.01 

42 

22 

* 

0.01 

23 

14 

* 

0.005 

49 

47 

* 

0.005 

29 

27 

* 

0.005 

46 

25 

* 

0.005 

27 

16 

* 

0.0025 

42 

47 

* 

0.0025 

24 

27 

* 

0.0025 

49 

28 

* 

0.0025 

29 

18 
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TABLE 6.10. Block diagonal preconditioner V (4.12) for FGMRES method 
(diagonal blocks are solved exactly) 


(a) Re 

= 1, Rm = 

1 

(b) Re = 

1, Rm = 

400 

(c) Re = 

400, 

Rm 

= i 

(d) Re = 400, 

Rm 

= 400 

h 

A 

!/4 

'A 

1/16 

k 

1/4 

'A 

1/16 

k 

At 

V 4 

Vs 

1/16 

k 

A 

1/4 

Vs 

!/ l 6 

0.02 

15 

16 

* 

0.02 

15 

16 

* 

0.02 

10 

6 

* 

0.02 

10 

6 

* 

0.01 

15 

14 

* 

0.01 

15 

13 

* 

0.01 

11 

7 

* 

0.01 

11 

7 

* 

0.005 

14 

13 

* 

0.005 

14 

13 

* 

0.005 

13 

8 

* 

0.005 

13 

8 

* 

0.0025 

12 

13 

* 

0.0025 

12 

13 

* 

0.0025 

15 

9 

* 

0.0025 

15 

9 
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Table 6.11. Block lower triangular preconditioner Me (4.23) for FGM¬ 
RES method (diagonal blocks are solved exactly) 


( a ) Re 

= 1 , 

Rm = 

1 

(b) Re = 

1 , R/?i = 

400 

( c ) Re = 

400, 

Rm 

= 1 

(d) Re = 400, 

Rm 

= 400 

k 

Af^\ 

1/4 
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1/16 

h 

1/4 

Vs 

1/16 


1/4 

Vs 
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/z 

1/4 

Vs 

1/16 

0.02 

56 

64 

89 

0.02 

36 

42 

70 

0.02 

39 

21 

14 

0.02 

24 

14 

10 

0.01 

54 

52 

74 

0.01 

34 

30 

49 

0.01 

43 

24 

17 

0.01 

25 

16 

12 

0.005 

50 

51 

56 

0.005 

29 

29 

38 

0.005 

47 

32 

21 

0.005 

27 

19 

14 

0.0025 

41 

50 

49 

0.0025 

25 

27 

26 

0.0025 

49 

38 

25 

0.0025 

32 

21 

16 


Table 6.12. Block diagonal preconditioner M (4.13) for FGMRES method 
(diagonal blocks are solved approximately) 


(a) Re 

= 1 , 

Rm = 

1 

(b) Re = 

1, Rm = 

400 

(c) Re = 

400, 

Rm 

= 1 

(d) Re = 400, 

Rm 

= 400 
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1/16 

h 
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1/4 
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1/16 
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18 
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11 

7 
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0.02 

11 

7 

5 

0.01 

16 

15 

21 
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16 

15 

20 

0.01 

13 

8 

6 

0.01 

13 

8 

6 

0.005 

16 

15 

17 
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16 

15 

17 
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15 

10 

7 

0.005 

15 

9 

7 

0.0025 

13 

15 

14 

0.0025 

13 

15 

14 

0.0025 

16 

11 

8 

0.0025 

16 

11 

7 


Table 6.13. Block lower triangular preconditioner Me (4.24) for FGM¬ 
RES method (diagonal blocks are solved approximately) 
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7. Conclusions 

In the numerical simulations for the incompressible MHD system, the most time- 
consuming part is solving the linear system. In this paper, we design several new block 
preconditioners for the linear systems resulting from the structure-preserving finite el¬ 
ement discretization proposed in [16]. We develop these preconditioners from the per¬ 
spective of the functional and PDE analysis. By rigorously proving the well-posedness of 
the discretization schemes and carefully studying the mapping property of the linearized 
operators between appropriate Sobolev spaces with proper weighted norms following 
the framework in [20, 19], we develop two types of preconditioners: norm-equivalent 
preconditioners (for example, block diagonal preconditioners) and FOV-equivalent pre¬ 
conditioners (for example, block triangular preconditioners). By revisiting the inf-sup 
conditions of the discrete systems under the new weighted norms, we theoretically ver¬ 
ify the robustness of the preconditioners when the time step size is small enough, and 
further prove that the resulting preconditioned Krylov iterative methods (for example, 
MINRES and GMRES methods) converge uniformly with respect to the mesh size, time 
step size, and physical parameters including the relative electrical conductivity <j r , the 
relative magnetic permeability }i r , the coupling number s, and the magnetic Reynolds 
number Rm. We also verify the theoretical conclusions by numerical experiments. 

Another contribution of this paper is that we improve the analysis of FOV-equivalent 
preconditioners proposed in [19]. Their analysis requires scaling parameters in front of 
the diagonal blocks in Me under certain constrains, which usually are difficult to choose 
in practice. In our analysis, with the help of an appropriate norm (•, •)» u we are able to 
remove those unnecessary scaling parameters, which makes the theoretical results con¬ 
sistent with practical implementation and observations. 

Finally we would like to point out that we have not considered the convection dom¬ 
inant cases in either the design of discretization or the design of preconditioner. For 
convection dominant cases (that could happen in (1.1) and or (1.3)), we need to use spe¬ 
cial discretization methods such as upwinding, streamline diffusion or other stabilization 
techniques and we also need to use tailored techniques for designing preconditioners for 
the resulting discrete systems. This is a topic of our ongoing investigations. We expect 
that for all practical purposes, preconditioners developed in this paper can be extended 
to certain convection dominant cases without much difficulty but we do not at all expect 
that it is easy to extend our theory to such cases as rigorous theoretical results that are 
uniform with respect to the magnitude of the convection are extremely rare for either 
continuous or discrete case in any circumstances. Despite the limitation of our results 
for the convection dominant case, we trust that the rigorous theoretical results and the 
supporting numerical experiments presented in this paper will help shed some light on 
certain aspects of the highly complicated MHD systems. 
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